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COMPUTATIONAL STUDY OF SHOCK WAVES PROPAGATING THROUGH 
AIR-PLASTIC-WATER INTERFACES. 

MAURICIO J. DEL RAZO* AND RANDALL J. LEVEQUE * 

Abstract. The following study is motivated by experimental studies in traumatic brain injury (TBI). Recent research 
has demonstrated that low intensity non-impact blast wave exposure frequently leads to mild traumatic brain injury (mTBI); 
however, the mechanisms connecting the blast waves and the mTBI remain unclear. Collaborators at the Seattle VA Hospital 
are doing experiments to understand how blast waves can produce mTBI. In order to gain insight that is hard to obtain by 
experimental means, we have developed conservative finite volume methods for interface-shock wave interaction to simulate 
these experiments. A ID model of their experimental setup has been implemented using Euler equations for compressible fiuids. 
These equations are coupled with a Tammann equation of state (EOS) that allows us to model compressible gas along with 
almost incompressible fiuids or elastic solids. A hybrid HLLC-exact Eulerian-Lagrangian Riemann solver for Tammann EOS 
with a jump in the parameters has been developed. The model has shown that if the plastic interface is very thin, it can be 
neglected. This result might be very helpful to model more complicated setups in higher dimensions. 
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1. Introduction. Recent research has demonstrated that low intensity non-impact blast wave exposure 
frequently leads to mild traumatic brain injury (mTBI) [11, 14]; however, the mechanisms connecting the 
blast waves and the mTBI remain unclear. Conventional imaging techniques, like computerized tomography 
and MRI, fail to reveal the damage when the injuries are mild. This suggests that the injury mechanisms 
might occur at very small length scales, even at the scale of a single cell. Some small scale damage could 
include metabolic cells’ changes or even alterations in the blood brain barrier(BBB) diffusion, as hinted in 
[ 8 ], 
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Fig. 1.1: (a) Cartoon of experimental setup. A cylindrical cross section of the shock tube with the sample is 
shown. The shock wave fronts travels through air in the shock tube hitting the plastic first, then the water 
and finally the endothelial cells sample, (b) Zoom-in of the Region R of (a). This region is the main concern 
of this work, and it will be modeled as one dimensional. 


In order to gain insight on the injury mechanisms, experimentalists at the Seattle Veterans Administra¬ 
tion hospital (VA hospital) are trying to reproduce the effects of blast waves in brain cells with the aid of a 
shock-tube. A sample with endothelial cells, the cells that form the BBB, is prepared and placed inside the 
shock-tube. The differences in several properties of the cells are studied before and after the blast wave has 
passed. The endothelial cells sample is placed inside of a polystyrene (plastic) cylindrical transwell container 
filled with an aqueous solution; a cartoon of the experimental setup is presented in Figure 1.1a. 
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The shock-tube is a long cavity filled with air, where the shock wave pressure profile can be measured. 
However, we don’t know what will happen when the shock wave interacts with the polystyrene cylindrical 
transwell filled with water. The physical effects of the shock wave crossing from air to plastic and from 
plastic to water are not evident from experimental data nor easy to obtain through experimental techniques. 
The model presented here concentrates on the one dimensional idealized scenario (Figure 1.1b) of the region 
R of Figure 1.1a. The computational simulations presented will help bring insight on the behavior of the 
shock wave crossing air-plastic-water interfaces, as the one shown in Figure 1.1b. They will also provide 
quantitative data of the strength of the shock wave before and after crossing the interfaces, and they will 
show that if the plastic interface is thin enough, it can be completely neglected. 

In Section 2, we show an exact result for linear acoustics of how the plastic interface becomes irrelevant 
as its thickness goes to zero. In Section 3, the appropriate model to represent the experimental setup will 
be presented; an overview of its numerical implementation is presented in Section 4. In Section 5, two 
computational experiments are described to show the plastic interface can be neglected without losing much 
accuracy. Moreover, we gained insights in the behavior of the shock wave when crossing the interfaces, 
of which the most relevant is the increase in pressure when the shock wave crosses from air to water. A 
discussion of the results is given on Section 6. These results form a basis for on-going work on this application, 
to be described further in future publications [2, 3]. 

2. Air-Plastic-Water interface for linear acoustics. Before trying to model the shock waves with 
complicated non-linear systems of conservation laws, we can start with a very simple case to gain insight: 
linear acoustics. In this case, we can compute the exact solution of the transmitted pressure through the 
air-plastic-water interface as a function of the acoustic impedance of each material and the plastic width. 
This can be derived from the fact that an acoustic wave with incident pressure po on the left of an interface 
between medium A (left) and H(right) produces a reflected and a transmitted wave with pressures given by 
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where Zk denotes the acoustic impedance of medium k. These relations can be easily derived from linear 
acoustics [9]. Now consider a triple one dimensional interface: air-plastic-water. With this setup, there 
will be an infinite number of reflections in the plastic layer. The transmitted wave in water after the 
reflection through the plastic yields 
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where Za^Zp^Z^ are the air, plastic and water impedances. Each transmitted wave increases the pressure 
behind the initially transmitted wave slightly and the asymptotic final amplitude of the transmitted wave is 
given by the sum of all these contributions, i.e: 
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Summing this geometric series yields 
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When the plastic layer is very thin this asymptotic value is quickly reached, and we note that it is exactly 
the same as if the plastic interface didn’t exist and we computed the transmission coefficient directly from 
air into water. Note we assumed the pressure profile on the left was a constant, po- However, this can be 
more complicated. It can have a decaying tail, in which case there will be interference from the tail in the 
reflected and transmitted waves. Nonetheless, assuming the plastic width is rco, the time elapsed between 
two transmitted waves in the water interface is given by r = 2n;o/cp, where Cp is the speed sound in plastic. 
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Therefore, as wq 0, the elapsed time does too. As a consequence, the interference from the tail will also 
disappear and the plastic interface can be neglected without losing accuracy. 

This calculation shows that for linear acoustics, if the plastic interface is very thin in comparison to 
the experiment characteristic length scales, the plastic interface can be neglected without significant loss of 
accuracy. From this result, we can expect similar results for the one-dimensional non-linear system that we 
will employ to model the experimental setup. 

3. The model. In order to perform a computational experiment, we first need to determine which are 
the most appropriate mathematical equations to describe the experimental setup. As we are interested in 
pressure waves (sound-waves), a first approach is to employ linear acoustic equations. However, as the form 
of the equations is linear, shock wave formation is not possible. A viable alternative is the Euler equations 
for compressible inviscid flow. As the equations are nonlinear, shock wave formation is accurately modeled. 
Furthermore, Euler equations also model conservation of energy, whose connection with temperature might 
be relevant to cell injury. Additionally, for this experimental setup, we are not concerned with large-scale 
movement of the fluid, so viscosity can be neglected by employing the inviscid equations. One more issue 
is that air has a high compressibility while water and plastic are almost incompressible. Using different 
parameters for the equations of state (EOS) for each material, we can model the three materials with the 
same equations. The equations are solved using the methods briefly explained in Section 4. 

A first approach to this problem is to use one-dimensional Euler equations by focusing in a cross section 
close to the central axis of Figure 1.1, 
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(3.1) 


where p is the density, u denote the velocities in the x (shock tube axis) direction, p is the pressure and E 
the internal energy. 

3.1. Tammann equations of state. The system of equations (3.1) is closed with the addition of an 
EOS. It is usually given as a relation between pressure, density and specific internal energy, i.e. p = p(p, e). 
The most well known EOS is the one for an ideal gas p = (7 — l)pe, where 7 is the heat capacity ratio, p 
is the density and e the specific internal energy. While this EOS is very good to describe the behavior of 
gasses, it might not be appropriate to model nearly incompressible materials like water or elastic solids. 

Several alternatives exist. In this work, we will use the stiffened gas EOS (SGEOS) or also known as 
Tammann EOS. This equation of state is very useful to model a wide range of fluids even in the presence of 
strong shock waves [5]. The Tammann EOS is given by 


p= ( 7 -l)pe-7 Poo, (3.2) 

where 7 and Poo can be determined experimentally for different materials. The Tammann EOS and the ideal 
gas EOS are the same except for the extra term — 7 Poo- As the pressure is a positive quantity, and 7 ,Poo > 0? 
this EOS can only be aimed for higher density fluids. For fluids with poo ^ 1, the relative change in density, 
when changing the pressure, is very small. Consequently, the Tammann EOS is a good approximation for 
fluids with a very low compressibility. It can also be helpful to model elastic solids, like plastic. It’s worth 
mentioning that for sufficiently weak shocks the Tammann EOS can be further simplified to the Tait EOS, 
see [5]. 

4. Numerical implementation. The Euler equations are a hyperbolic system of conservation laws, so 
they can be solved employing finite volume methods (FVM). This is done by employing the wave propagation 
algorithms described in [9] and implemented in Clawpack [ 1 ]. The fundamental problem to solve at each cell 
interface of our computation is the well known Riemann problem. The equations of motion are solved by 
implementing a hybrid Riemann HLLC-exact type approximate solver for one-dimensional Euler equations 
with interfaces. This solver couples a Eulerian HLLC approximate Riemann solver to an exact Riemann 
solver for the Tammannn EOS in Lagrangian coordinates. As the interfaces are represented by contact 
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discontinuities, the HLLC solver is ideal to deal accurately with interface problems. The method can easily 
be extended to two dimensions by employing dimensional splitting. An overview of the HLLC Riemann 
solver and the coupling is given next. 

4.1. The HLLC solver. A general one-dimensional Riemann problem for a system of conservation 
laws like Euler equations can be stated as 

Q + F(Q), = 0. Q(*.0) = { (4.1) 

where Q is the vector of state variables. The HLLC (Harten-Lax-van Leer-Contact) solver is an approximate 
Riemann solver to the problem (4.1). As the Riemann solution to the one-dimensional Euler equations 
consists of three waves: two acoustic waves with a contact discontinuity in between. The main idea of the 
HLLC solver is, given upper bounds for the left and right going wave speeds Sl and Sr^ assume a wave 
configuration of three waves separating four constant states. The approximate solution for this method will 
be of the form 


Q(x,t) 


^ Ql 
Q*l 

, Qi? 


if f < 

if 5 'l < f < S'*, 
if S* < f < Si,, 
if f > Si,, 


where is the approximate wave speed of the contact discontinuity. Assuming we can obtain Sr and Sr by 
some other algorithm, we only need to find Q*l, Q*i? and S* to solve the problem. These quantities can be 
obtained by integrating over a box in the x, t plane, using the Rankine-Hugoniot conditions and assuming 
constant pressure and normal velocity across the contact discontinuity, see [13]. The desired states for the 
Euler equations with Q = [p^pu^EY are given by 


Q*/c 


S/cQ/e F 
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with: D = [0,1, S*]^, 


(4.2) 


where the contact discontinuity speed is given by 


Pr — Pl PlUl{Sl — ul) — PrUr{Sr — ur) 
Pl{Sl — Ul) — Pr{Sr — Ur) 


and the pressure in the middle state p* is 


P^ = Pl E Pl{Sl — ul){S^ — ur) = PrE Pr{Sr — ur){S^ — ur) 

where p/^, with k = L^R are the left or right density and speed in the Euler equations, for detailed 
calculations the reader is referred to [13]. 

In order to calculate the wave speeds Sr and Sr, we will need to calculate the sound speed. This is 
where we require the EOS. A simple estimate is the one given by Davis [13] as Sr = min{rtx. “ ur — cr} 
and Sr = max{rti, E cr,ur cr}, where Uk is the normal velocity and Ck is the sound speed on each side, 

k = L,R. The speed of sound using the Tammann EOS is given by c = A possible improvement 

is to employ Roe averages in wave speed estimates [4] . These Roe averages can be calculated using different 
configurations, some might be more accurate when dealing with interfaces, as pointed out in [6] 

The HLLC solver just discussed works well for solving the one dimensional Euler equations with an ideal 
gas EOS. However, we want to implement the HLLC solver with the Tammann EOS across an air-plastic- 
water interface. The difference between the parameters for different materials in the Tammann EOS is of 
several orders of magnitude, as shown in Table 5.1. This generates instabilities in the HLLC solver, more so 
in a two or higher dimensional setting. The instability is due to the fact that we model the interfaces as being 
fixed in space; however, there is always a displacement of the contact discontinuity, i.e. the interface. In order 
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to solve this issue, we model the air in Eulerian coordinates using the usual HLLC solver. If any of the cells 
is a water or plastic cell, we modify our original HLLC or exact solver to work in Lagrangian coordinates, 
where the interface is actually fixed with respect to the reference frame. This is done by displacing the frame 
of reference by S'*, 


Sl = Sl-S, S* = 0 Sr = Sr-S, (4.3) 

For instance, assume we are running a one dimensional simulation of the Euler equations, with a fixed 
interface modeled by a jump in the parameters of the EOS. The interface is aligned to the edge between 
cells i and i +1, the transformed Riemann solver will be as shown in Figure 4.1. This will ensure the contact 
discontinuity velocity is zero and consequently the interface is modeled as fixed. The wave contributions 
will be the correct ones, since we are just modifying the wave velocity and not the solution g’s. There is, of 
course, an error made at the interface when coupling the two descriptions; however, as the displacements of 
the interface are very small due to very low compressibility, this error is not big enough to cause instabilities 
as before. The speeds (4.3) along with the solution in the different regions (4.2) are the output of our 
Riemann solver into the Finite Volume method implemented into Clawpack [1]. Note the transformed 
speeds in (4.3) should only be used in the presence of an interface; otherwise, the Riemann solver should 
output the non-transformed speeds and Sr. 




Fig. 4.1: Transformation from Eulerian coordinates to Lagrangian coordinates for the HLLC or exact Rie¬ 
mann solver between grid cells i and i + 1. The transformation is only done when one or more of the cells 
are part of the highly incompressible material where S^ <^1. 


In order to provide better accuracy in certain regions, we also developed and implemented an exact 
Riemann solver for the Euler equations using the Tammann EOS with a jump in the parameters, similar to 
the one by Ivings & Toro [7]. For the sake of brevity, this will be presented in a forthcoming manuscript. 
The HLLC solver will be used to model the air in Eulerian coordinates, and the exact solver will be used to 
model the interface and the almost incompressible material in Lagrangian coordinates. The transformation 
to Lagrangian coordinates for the exact solver is exactly the same one as in equations (4.3). 

5. Computational experiments. We want to study the importance of the plastic interface in the 
computational experiment illustrated in Figure 1.1. The simplest scenario is a one-dimensional model of the 
Region R of the transwell in Figure 1.1b, i.e an air-plastic-water interface. In this section, we will study a 
shock wave crossing a general air-plastic-water interface for different plastic widths. This will bring insight 
on the behavior of the shock wave, and it will show how relevant is the thin plastic interface. We will 
analyze the one-dimensional air-plastic-water interface problem, and then we’ll compare it to the simpler 
one-dimensional air-water interface. 

5.1. Air-plastic-water interface. The first step is to input the right initial conditions into our sim¬ 
ulation. The actual form of the shock wave traveling through the shock tube before hitting the transwell 
was obtained experimentally. The sensor outputs pressure amplitude as a function of time. Assuming an 
average speed of sound in air, it can be converted to a function of distance as shown in Figure 5.2. The shape 
can be broadly approximated by an idealized shock wave (dashed line in Figure 5.2). This approximated 
form of the shock wave is introduced as part of the initial condition on the left part of the domain in the 
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(a) 


(b) 


(c) 


Fig. 5.1: Shock wave crossing the air-plastic-water interface at different times. The arrows indicate the 
position of the 4 gauges that measure the pressure as a function of time. The gauges are numbered from left 
to right, and the plastic interface width for this case is 2 . 6 m 


simulation; however, this is not a trivial procedure since we must input the density, momentum and energy, 
and we only have the pressure. Using the isentropic EOS, the ideal gas EOS and the expression for the 
speed of sound, an educated guess for the initial condition in terms of the pressure is given far away from the 
transwell, elsewhere the initial condition will be constant: constant density, zero velocity and atmospheric 
pressure. This initial condition is then modified until the right amplitude and shape of the shock wave front 
is obtained. The resulting shape of the shock wave before hitting the interface can be seen in Eigure 5.1a. 
The pressure is measured in KPa with an ambient base pressure of lATM = 101.325KPa. 

The one-dimensional equations with the triple interface are solved using the methods mentioned in 
Section 3. The different materials are modeled using different parameters for the Tammann EOS, see Section 
3.1. The choice of parameters is shown in Table 5.1. 


Material 

7 

Poo {G Pa) 

Air (Ideal gas EOS) 

1.4 

0.0 

Plastic (polystyrene) 

1.1 

4.79 

Water 

7.15 

0.3 


Table 5.1: Parameters for the Tammann EOS to model the different materials. The parameters for air and 
water were taken from [5]. Since the polystyrene is a solid, 7 was chosen to be close to 1 , and was 
adjusted to yield the right speed of sound in polystyrene. 


The solution of the shock wave crossing the three interfaces; air, plastic and water, at different times 
is shown on Eigure 5.1. It can be observed that every time the pressure wave hits an interface, part of the 
wave is reflected and part of it is transmitted. This effect can occur multiple times depending on how the 
interfaces are set up. We can also observe that the amplitude of the shock wave increases as it passes from 
air to plastic and decreases when passing from plastic to water. This effect is mostly due to the continuity 
of pressure at the change in compressibility. In order to keep the pressure at the interface continuous, the 
transmitted wave amplitude has to be the same as the sum of the incident wave and the reflected wave. When 
the compressibility is very high in the adjacent material, the interface will behave similarly to a solid wall. In 
this case, since the reflected wave will have an amplitude almost equal to the incident wave, the transmitted 
wave could have an amplitude almost twice as big as that of the incident wave. This explains why the 
pressure can increase or decrease when crossing an interface. As a consequence of the different behaviors in 
the interface, even for the one-dimensional case, a complex behavior is observed. These numerical simulations 
provide accurate insight in situations where simple intuition might be insufficient. 
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Width (m) 

luitial(KPa) 

Gauge 2 (KPa) 

Gauge 3 (KPa) 

Gauge 4 (KPa) 

2.6 

184.06 

247.76 

305.88 

258.24 

1.4 

184.06 

207.53 

298.19 

259.71 

0.6 

184.06 

187.72 

283.72 

274.90 

0.2 

184.06 

183.31 

282.34 

280.18 

0.1 

184.06 

183.31 

284.29 

283.55 


Table 5.2: The maximum amplitude measured at three pressure gauges for differeut widths of the plastic 
iuterface. The iuitial shock wave is the same for all cases, aud the gauge plots are placed before, iuside aud 
after the plastic iuterface as showu iu Figure 5.1. 



Fig. 5.2: The shock wave form obtaiued from the seusor iuside the shock tube is showu as the solid thiu hue. 
The coarse approximatiou to be used as au iuitial couditiou iu our simulatiou is showu with a dashed hue. 
Au average speed of souud of c = 344 m/s is assumed. 


Iu Figure 5.2, we show from experimeutal data the shock wave profile iu the air before hittiug the 
trauswell; however, we are iuterested iu the shape aud amplitude of the shock wave iu the water just before 
hittiug the eudothelial cells sample. Iu order to do so, we first ueed to kuow how relevaut the plastic iuterface 
is iu our model, or if it is relevaut at all. Computatioually, the plastic iuterface is hard to model because the 
width of the plastic is very small ( mm) iu comparisou to the shock wave waveleugth (meters). The followiug 
experimeut explores how the width of the plastic iuterface affects the shock wave profile. Additioually, we 
show au accurate model cau be obtaiued eveu wheu completely iguoriug the plastic iuterface. 

Iu Figure 5.1, a set of four gauges is showu that allow to measure the pressure at certaiu poiuts as a 
fuuctiou of time. The maximum amplitude of the pressure profile was measured at gauges 2,3 aud 4 for 
differeut widths of the plastic iuterface. The plastic is always assumed to be ceutered at x = 0. The results 
are preseuted iu Table 5.2. Iu Figure 5.3, the full pressure profiles as a fuuctiou of time are showu at the 
three gauges for three of the plastic widths showu iu Table 5.2. 

The results iu Table 5.2 aud Figure 5.3 show the maximum amplitude at gauge 2 is reduced as the plastic 
width is decreased. Not surprisiugly, this is a cousequeuce of haviug less iuterfereuce with the reflected shock 
wave, siuce the gauge is farther away from the iuterface as the plastic width is reduced. This effect is clearly 
showu ou Figures 5.3a, 5.3d, 5.3g. The maximum amplitude at gauge 3 is somewhat dimiuished at first; 
however, it seems to be reachiug a plateau arouud 280.0KPa. The behavior at gauge 3 is uot trivial; the 
shock wave bouuces back aud forth several times, iuterferiug with itself coustautly. Iu Figures 5.3b, 5.3e, 
5.3h, we cau see the iuterfereuce becomes so fast that the pressure profile iu the plastic seems to couverge 
to a shock wave shape as the plastic width is reduced. At gauge 4, we cau observe the iuterfereuce betweeu 
the set of trausmitted shock waves geuerated by the back aud forth reflectious withiu the plastic iuterface. 
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As the plastic width is reduced, the time elapsed between the transmitted shock waves is reduced and the 
interference increased. Nonetheless, when the plastic width is very small, the interference becomes so fast 
that the pressure profile seems to converge again to a shock wave shape, as shown in Figures 5.3c, 5.3f, 
5.3i. Furthermore, note the difference in the shock wave shape in Figures 5.3h, 5.3i is almost unnoticeable. 
It almost seems like the shock wave is only crossing one interface instead of two. This motivates the next 
experiment. 



(a) (b) (c) 





(d) 


(e) 


(f) 





(g) 


(h) 


(i) 


Fig. 5.3: Pressure (KPa) gauge plots as a function of time (seconds). Each row of figures shows the three 
gauge plots for three different widths (2.6m, 1.4m and 0.1m) of the plastic interface, as shown in Table 5.2. 
The plots (g) and (i) for gauge 2 and 4 also show the pressure gauge plots when there is no plastic interface 
at all; the difference is almost unnoticeable. 


5.2. Air-water interface. In reality, the plastic is so thin that is really unnoticeable on larger scales. 
Furthermore, as the plastic is almost an incompressible medium, one should expect it would transfer the 
shock wave infinitely fast without energy loss. Therefore, instead of the triple material interface, let’s now 
consider only an air-water interface. The results of this simulation are shown in Figure 5.4. The gauge plots 
for gauge 2 and 4 are shown in Figures 5.3g and 5.3i, along with the thin plastic results. The maximum 
amplitude in each of these gauges is presented in Table 5.3. 
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(a) (b) 

Fig. 5.4: Shockwave before and after crossing the air-water interface. The arrows indicate the position of the 
gauges that measure the pressure as a function of time. The gauges are the same as in Figure 5.1; however, 
gauge 3 is not shown for plotting aesthetics. 


Initial (KPa) 

Gauge 2 (KPa) 

Gauge 4 (KPa) 

184.06 

184.40 

284.26 


Table 5.3: The maximum amplitude measured by two pressure gauges before and after the air-water interface. 
The initial shock wave is the same as in the air-plastic-water case. 


Comparing the air-water interface results against the ones for the smallest plastic width in the air-plastic- 
water interface case, we can observe the percentage error in the maximum pressure amplitude of gauge 4 is 
of 0.38%. This is also obvious from the thin plastic and no plastic comparison in Figures (5.3g, 5.3i). An 
analytic result of this effect for linear acoustics is shown in Section 2. These results will allow us to simplify 
higher dimensional air-plastic-water interface problem to a simpler air-water interface. However, it should 
be noted that these results for one dimensional interfaces might not be as accurate when employed in higher 
dimensions. 

6. Discussion. Using computational experiments, we showed that there is not a significant difference in 
the transmitted shock wave between the air-plastic-water and the air-water interface. Therefore, thin plastic 
interfaces can be neglected in future computations. We also observed an amplification and elongation of the 
shock wave when crossing the interface. This effect is due to the different sound speed and compressibility 
of the materials. The amplitude of the initial pressure wave in air increased by 54% when measured in 
the water. Since the shock intensity is much higher in water, where the endothelial cells sample is, the 
amplification effect is highly relevant to study injury mechanisms. Finally, the amplification effect also 
occurs when passing from air to plastic or any solid material, like bone. This means that a pressure shock 
wave perceived as weak outside in air might be much more intense when traveling through the skull and the 
brain. Future work will involve an extension to three dimensions and implementation of two-phase models 
like those of Pelanti & Saurel [10, 12]. 
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